CAM Implementation of Mesoscale Eddy Diffusivity in Terms of Mixing Length Theory
MICHAEL BATES, ROSS TULLOCH, JOHN MARSHALL, AND RAFFAELE FERRARI
Steering Level status
1. L_eddy values ok
2. Rossby wave speed ok
3. Eady inverse timescale
a.has anomalous equatorial belt (f goes to 0 at equator)
b.Richardson averages are too high causing sigma to drop too low
4. Zonal Eddy Phase Speed (c)
a. Too high 10S:10N
b. Missing westward velocities in ACC which appear
in Hughes Data used in this parameterization
5. (U-c) Because of problem with c this term is also too high 10S:10N.
6. U_rms (alpha*sigma*l_eddy)
a. Overall too weak (due to problems with Eady timescale.)
b. Missing alot of structure especially in the NH. Western Boundary currents missing
some wbc structure in eady but multiplying by leddy removes that at higher latitudes
7. Zonal mean velocity looks ok compared to ECCO annual average
8. Suppression, Lmix, and K are wrong because of errors in the above terms
Todo
1. New calculation for sigma may correct sigma deficiencies and fix the U_rms problems.
or fixing equator with max (/f/,sqrt(c_r*2*Beta)) and changing calculation for sqrt(R_i) would improve sigma
2. Need to work on Zonal phase speed c
3. Fill in rest of plots for K
4. Plot vertical distribution of K
%pylab inline
from IPython.display import display
from IPython.display import Image
from netCDF4 import Dataset
from netCDF4 import num2date
import numpy as np
import sys
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from matplotlib.colors import LogNorm
from matplotlib.ticker import LogFormatter
from mpl_toolkits.basemap import Basemap, shiftgrid, addcyclic
#for interpolation (you will have to install pyresample first)
matplotlib.rc('text', usetex=True)
matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
import pyresample
numpy.version.version
def make_basemap(projection='cyl', resolution='c',lon_0=0, ax=None,llcrnrlon=25):
m = Basemap(projection=projection, resolution=resolution,
ax=ax,llcrnrlat=-90,urcrnrlat=90,
llcrnrlon=llcrnrlon,urcrnrlon=llcrnrlon+360,fix_aspect=False)
m.drawcoastlines()
m.drawmapboundary(fill_color='white')
m.fillcontinents(color='k',lake_color='white')
parallels = np.arange(-60, 90, 30.)
meridians = np.arange(-360, 360, 60.)
m.drawparallels(parallels, labels=[1, 0, 0, 0])
m.drawmeridians(meridians, labels=[0, 0, 1, 0])
return m
def MapZoneContour(lon,lat,data,figsize=(14,9),title=None,levels=None,fmt=None,norm=None):
map2zonal=6
fig = plt.figure(figsize=figsize)
ax1 = plt.subplot2grid((1,map2zonal), (0,0), colspan=(map2zonal-1))
m=make_basemap()
#make longitudes monotonic for contouring add longitude period to end
lon_mono = np.where(np.greater_equal(lon,min(lon[:,0])),lon-360,lon)+360.
endlon=lon_mono.shape[1]-1
lon_mono[:,endlon]=lon_mono[:,endlon]+360.
# contour and labels
if (levels==None):
CS=m.contour(lon_mono,lat_targcyc,data,10, linewidth=.5, colors='k')
else:
CS=m.contour(lon_mono,lat_targcyc,data,levels=levels, linewidth=.5, colors='k')
plt.clabel(CS, CS.levels,fmt=fmt)
plt.suptitle(title,fontsize = 15,y=1.04)
# zonal subplot
ax2 = plt.subplot2grid((1,map2zonal), (0,(map2zonal-1)), sharey=ax1)
# First find the zonal mean by averaging along the latitude circles
data_zonal = data.mean(axis=1)
plt.plot(data_zonal,lat[:,0])
return fig
def MapContourf(lon,lat,data,
figsize=(14,9),
title=None,
levels=None,
fmt=None,
extend='neither',
contourlines=False,
addzonal=False,
invertcb=False,
norm=None,
yscale='Linear'):
# from matplotlib.colors import LogNorm
map2zonal=6
formatter = ticker.ScalarFormatter()
formatter.set_scientific(True)
formatter.set_powerlimits((0,0))
fig = plt.figure(figsize=figsize)
#make longitudes monotonic for contouring add longitude period to end
lon_mono = np.where(np.greater_equal(lon,min(lon[:,0])),lon-360,lon)+360.
endlon=lon_mono.shape[1]-1
lon_mono[:,endlon]=lon_mono[:,endlon]+360.
if (fmt==None):
fmt='%.3e'
if (addzonal):
ax1 = plt.subplot2grid((1,map2zonal), (0,0), colspan=(map2zonal-1))
else:
ax1 = fig.add_subplot(111)
if (contourlines):
# Create contour lines
CS0 = ax1.contour(lon_mono,lat,data,levels=[0.],norm=norm,linewidths=0.5,colors='k',latlon=True)
# y-labels on the right
ax1.yaxis.tick_right()
if (invertcb):
CS1 = ax1.contourf(lon_mono,lat,data,levels=levels,norm=norm,cmap='jet_r')
cbar = fig.colorbar(CS1, ax=ax1,format=fmt)
cbar.ax.invert_yaxis()
if (extend!=None):
CS1.cmap.set_under('darkred')
CS1.cmap.set_over('black')
else:
if (norm==None):
CS1 = ax1.contourf(lon_mono,lat,data,levels=levels,extend=extend,norm=norm,cmap='jet')
cbar = fig.colorbar(CS1, ax=ax1, format=fmt)
else:
CS1 = ax1.contourf(lon_mono,lat,data,levels=levels,norm=norm,cmap='jet')
l_f = LogFormatter(10, labelOnlyBase=False)
cbar = fig.colorbar(CS1, ax=ax1, ticks=levels, format=l_f)
if (extend!=None):
CS1.cmap.set_under('black')
CS1.cmap.set_over('darkred')
m=make_basemap()
plt.suptitle(title,fontsize = 15,y=1.04)
fmt = matplotlib.ticker.FormatStrFormatter("%.3f")
if (addzonal):
# zonal subplot
ax2 = plt.subplot2grid((1,map2zonal), (0,(map2zonal-1)), sharey=ax1)
formatter1 = ticker.ScalarFormatter()
formatter1.set_scientific(True)
formatter1.set_powerlimits((0,0))
ax2.xaxis.set_major_formatter(formatter1)
# First find the zonal mean by averaging along the latitude circles
data_zonal = data.mean(axis=1)
if (levels!=None):
pylab.xlim([levels.min(),levels.max()])
plt.plot(data_zonal,lat[:,0])
return fig
basedir='/scratch/jet'
### get Variables from file
#nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.004.pop.h.nday1.0249-01-02.nc')
#nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.006.pop.h.nyear1.0249.nc')
nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.yr.0249.nc')
#print(nc.variables)
sigma_avg = nc.variables['SIGMA_AVG'][0,:,:]
utgrd_max = nc.variables['UTGRD_MAX'][0,:,:,:]
vtgrd_max = nc.variables['VTGRD_MAX'][0,:,:,:]
stgrd_max = nc.variables['STGRD_MAX'][0,:,:,:]
ttgrd_max = nc.variables['TTGRD_MAX'][0,:,:,:]
###sigma_gamma = nc.variables['SIGMA_GAMMA'][0,:,:]
###maxftinv = nc.variables['MAXFTINV'][0,:,:]
fcor = nc.variables['FCOR'][0,:,:]
c_eddy = nc.variables['C_EDDY'][0,:,:]
c_rossby = nc.variables['C_ROSSBY'][0,:,:]
l_eddy = nc.variables['L_EDDY'][0,:,:]
l_mix = nc.variables['LMIX'][0,:,:]
l_rhines = nc.variables['L_RHINES'][0,:,:]
l_rossby = nc.variables['L_ROSSBY'][0,:,:]
l_rossbyeq = nc.variables['L_ROSSBYEQ'][0,:,:]
urms_avg = nc.variables['URMS_AVG'][0,:,:]
urmsavgsq = nc.variables['URMS_AVG_SQ'][0,:,:]
uminuscsq = nc.variables['U_MINUS_C_SQ'][0,:,:]
suppress = nc.variables['SUPPRESS'][0,:,:]
umean = nc.variables['U_MEAN'][0,:,:]
umean10 = nc.variables['U_MEAN'][0,:,:]
#umean50 = nc.variables['U_MEAN_50'][0,:,:]
#umean100 = nc.variables['U_MEAN_100'][0,:,:]
#umean200 = nc.variables['U_MEAN_200'][0,:,:]
#umean500 = nc.variables['U_MEAN_500'][0,:,:]
kappalat = nc.variables['KAPPA_LATERAL'][0,:,:]
#tinv = nc.variables['TINV'][0,:,:]
grate = nc.variables['GRATE'][0,:,:]
btp = nc.variables['BTP'][0,:,:]
hmxl2 = nc.variables['HMXL_2'][0,:,:]
xmxl2 = nc.variables['XMXL_2'][0,:,:]
lat_orig = nc.variables['TLAT'][:]
lon_orig = nc.variables['TLONG'][:]
###time = nc.variables['time']
dz = nc.variables['dz'][:]
dzw = nc.variables['dzw'][:]
nc.close()
nc_ecco = Dataset(basedir+'/steering/ecco_ann_zonal.nc')
ecco_u=nc_ecco.variables['U'][0,0,:,:]
lat_ecco = nc_ecco.variables['LATITUDE_T'][:]
lon_ecco = nc_ecco.variables['LONGITUDE_UN179_181'][:]
nc_ecco.close()
###time_conv = num2date(time, time.units)
# this is the lefthand longitude starting point for plotting
llcrnrlon=25
# proportion of plot to side zonal average ie. 6:1
map2zonal=6
### use time_conv for plots using time as an axis
###fig = figure()
###plot(time_conv, air[:,10,10])
###fig.autofmt_xdate()
### or ###
###contourf(lat[:],time_conv,air[:,:,10])
###colorbar()
# setup for using pyresample to interpolate plot
# make sure lon is in range of -180:180 for pyregrid
lon_orig = (lon_orig + 180) % (360) - 180
# open dataset for output grid get latout,lonout
fc = Dataset(basedir+'/steering/co2flux_fossil_1751-2006-monthly_0.9x1.25_c20100204.nc')
lon= fc.variables['lon'][:]
lat = fc.variables['lat'][:]
fc.close()
# get 2D versions of the lat and lon variables add longitude start point here!
lon_targ, lat_targ = np.meshgrid(lon+llcrnrlon, -1.*lat)
# add cyclic point to output lat,lon arrays for plot routines
lon_targcyc = np.concatenate((lon_targ, lon_targ[:, 0:1]), axis=1)
lat_targcyc = np.concatenate((lat_targ, lat_targ[:, 0:1]), axis=1)
# make sure lon is in range of -180:180 for pyregrid
lon_targcyc = (lon_targcyc + 180) % (360) - 180
#When all the necessary inputs are read in memory, that's where we start using pyresample's magic.
#Create a pyresample object holding the origin (POP) grid:
orig_def = pyresample.geometry.GridDefinition(lons=lon_orig, lats=lat_orig)
#Create another pyresample object for the target (curvilinear) grid:
targ_def = pyresample.geometry.GridDefinition(lons=lon_targcyc, lats=lat_targcyc)
print len(lon_targcyc[0,:])
print len(lat_targcyc[:,0])
omega=7.2921e-5
deg2rad=pi/180
f=2*omega*sin(lat*deg2rad)
Update the matplotlib configuration parameters:
mpl.rcParams.update({'font.size': 18, 'font.family': 'serif'})
Implementation Notes on Diffusivity with Steering Level Suppression
12 November 2014
Starting with
$$K=u_{rms}∗L_{mix}$$
the general form for diffusivity, $K$, is given by equation (6) of Bates et al. (2014). Namely,
$$L_{mix} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $$
Here,
$u_{rms}$ is the root-mean-square (rms) eddy velocity;
$L_{mix}$ is the mixing length;
$L_{eddy}$ is the eddy diameter (depth independent);
$u_{mean}$ is the mean zonal velocity (resolved);
$c$ is the zonal eddy phase speed (depth independent);
$\Gamma = 0.35$;
$b1 \sim 4$.
Note that
$L_{mix} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}$
Eddy Length Scales:
Rossby deformation radius = $L_r = {c_r \over |f|}$
Equatorial Rossby deformation radius $= L_{req} = {\sqrt{c_r \over 2\beta}}$
Rhines scale $= L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}$
Where, $c_r$ is the first baroclinic wave speed computed following equation (2.2) of Chelton et al. (1998) with $m=1$;
$f$ is the Coriolis parameter;
$\beta$ is the latitudinal variation of the Coriolis parameter; and
$\sigma_{vi}$ is the Eady growth rate given by
$\sigma_{vi} = {f \over \sqrt{R_i}}$
with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.
So, any one these length scales could be used as an eddy length scale. An alternative is
$L_{eddy} = min (L_r, L_{req}, L_{Rh}).$
Eddy Velocity:
$u_{rms} = alpha*\sigma_{vi}*L_r$
where $\sigma_{vi}$ is the Eady growth rate based on local Richardson number and $\alpha$ is a scaling constant.
Zonal phase speed:
$c = - \beta * L_r^2$
Bates, M., R. Tulloch, J. Marshall, and R. Ferrari, 2014: Rationalizing the spatial distribution
of mesoscale eddy diffusivity in terms of mixing length theory. J. Phys. Oceanogr., 44,
1523-1540, doi: 10.1175/JPO-D-13-0130.1.
Chelton, D. B., R. A. deSzoeke, M. G. Schlax, K. E. Naggar, and N. Siwertz, 1998: Geographical
variability of the first baroclinic Rossby radius of deformation. J. Phys. Oceanogr., 28,
433-460.
Tullock, R., J. Marshall, and K. S. Smith, 2009: Interpretation of the propagation of surface
altimetric observations in terms of planetary waves and geostropic turbulence. J.
Geophys. Res., 114, C02005, doi: 10.1029/2008JC005055.
Questions:
In the 2D implementation, $u_{rms}$ and $u_{mean}$ specifications: upper-ocean vertically or integrated or at $z = 0$?
A: $U_{rms}$ is not depth dependent; both $u_{rms}$ and $u_{mean}$ are for surface only
In the 2D implementation, vertical profile will be specified by $N2(z)$?
A: Yes, $N^2(z) \over N_{ref}(z)$ to be more precise
Local $R_i$ use imbedded in sigma in $u_{rms}$ calculation?
A: No, $u_{rms}$ is depth independent
alpha = ?
A: trial and error
Cancellation of $f$’s in $u_{rms}$ calculation?
Zonal phase speed equation correct? Both $\beta$ and $L_r$ will be positive, producing $c < 0$ always. This appears to be in contrast with Tullock et al. (2009).
A: What is plotted in Bates is (U-c)
An alternate derivation of $\sigma_{vi}$ is:
$ R_i = {N^2 \over { ( \frac {\partial u}{\partial z} )^2 + ( \frac {\partial v}{\partial z} )^2} }$
$ N^2 = {{-g \over \rho_0 }\frac {\partial \rho}{\partial z}} $
After hydrostatic and geostrophic approximations
$f \frac {\partial v}{\partial z} = {{-g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad f \frac {\partial u}{\partial z} = {{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $
so
$\frac {\partial v}{\partial z} = {{-1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad \frac {\partial u}{\partial z} = {{1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $
$\therefore$
$ R_i = {f^2 N^2 \over \underbrace{{ {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial y})^2 + {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial z})^2} }_\text{m^4}}\quad\quad = \quad\quad {f^2N^2 \over m^4} $
$\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}$
$RX_1 = RX_{east} = \Delta\rho_x = \rho_{i+1,j} - \rho_{i,j}$
$RY_1 = RY_{north} = \Delta\rho_y = \rho_{i,j+1} - \rho_{i,j}$
$RZ_1 = RZ_{k+1} = \Delta\rho_z = \rho_{k} - \rho_{k+1}$
$\displaystyle{1 \over L_{R_i}} \displaystyle\int_{2000m}^{100m} \left\lbrace { {-d\over\rho_0}{\frac {\partial \rho} {\partial z}} \over { {g^2 \over \rho_0^2 } \left[( \frac{\partial \rho}{\partial y})^2 + ( \frac{\partial \rho}{\partial z})^2\right] } } \right\rbrace dz$
Note: missing $f^2$ which will be cancelled when forming $\sigma_{vi}$
$\quad\quad$ so $\cdots$ this is not $R_i$
Implementation notes
For Level K
Numerator : Top $= -grav * RZ_{SAVE}(\cdots k+1) * dzwr(k)$
Denominator :
$ \begin{align} work1 = p25 & * ( RX(..,i_{east},k)^2 \\ & + RX(..,i_{west},k)^2 \\ & + RX(..,i_{east},k+1)^2 \\ & + RX(..,i_{west},k+1)^2 ) / DXT(i,j)^2 \\ \end{align} $
$ \begin{align} work2 = p25 & * ( RY(..,j_{north},k)^2 \\ & + RY(..,j_{south},k)^2 \\ & + RY(..,j_{north},k+1)^2 \\ & + RY(..,j_{south},k+1)^2 \\ \end{align} $
$work3 = {\left( TOP \over (grav^2*(work1+work2))\right)}*dzw(k)$
Notes:
1)Need to be careful at top and bottom of ocean
2)Accurate dzw(k) for each (i,j) to form $L_{R_i}$
3) When constructing $sigma$ itself, use $RZ_{SAVE}$ with a minimum N value
4) use eps2
CESM First Baroclinic Wave Speed $c_r \cancel {= {1\over\pi} \int_{-2000}^{-100} N(z)dz}$ Now it is over full model depth
#Resample (aka re-project, re-grid) the NCEP data to target grid. First with nearest neighbour resampling...
c_rossby_nearest = pyresample.kd_tree.resample_nearest(orig_def, c_rossby, \
targ_def, radius_of_influence=500000, fill_value=None)
fig1=MapZoneContour(lon_targcyc,lat_targcyc,c_rossby_nearest,figsize=(14,9),
title="CESM First Baroclinic Wave Speed (cm/s)",
fmt='%.1i' )
c_r_chelton=Image(filename=basedir+'/steering/chelton_sfig1.jpg')
display(c_r_chelton)
First Baroclinic Rossby Radius Chelton vs CESM
CESM Rossby deformation radius = $L_{eddy}$
$L_r = {c_r \over |f|}$
$L_{req} = {\sqrt{c_r \over 2\beta}}$
$\require{cancel} \cancel {L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma \over \beta}}$
$L_{eddy} = min (L_r, L_{req}, \cancel{L_{Rh}}).$
#####Resample (aka re-project, re-grid) the NCEP data to target grid. First with nearest neighbour resampling...
####l_eddy_nearest = pyresample.kd_tree.resample_nearest(orig_def, l_eddy, \
#### targ_def, radius_of_influence=500000, fill_value=None)
fcor_nearest = pyresample.kd_tree.resample_nearest(orig_def, fcor, \
targ_def, radius_of_influence=500000, fill_value=None)
absf=np.abs(fcor_nearest)
print np.max(absf),np.min(absf)
l_rossby_calc=c_rossby_nearest/absf
print l_rossby_calc.max(),l_rossby_calc.min()
btp_nearest = pyresample.kd_tree.resample_nearest(orig_def, btp, \
targ_def, radius_of_influence=500000, fill_value=None)
l_rossbyeq_calc=np.sqrt(c_rossby_nearest/(2*btp_nearest))
l_eddy_calc=np.minimum(l_rossby_calc,l_rossbyeq_calc)
### convert regridded radius from cm to km
cm2km=1.e-5
print np.min(l_eddy_calc)
l_eddy_calc_km=l_eddy_calc*cm2km
clevs=[10,20,30,40,50,60,80,100,150,230]
fig2=MapZoneContour(lon_targcyc,lat_targcyc,l_eddy_calc_km,figsize=(14,9),
levels=clevs,
title="CESM First Baroclinic Rossby Radius (km) $L_{eddy} \sim min(l_r,l_{req})$",
fmt='%.1i' )
l_r_chelton=Image(filename=basedir+'/steering/chelton_fig2.jpg')
display(l_r_chelton)
$\sigma$ Eady growth rate
$\sigma$ is the Eady growth rate given by
$\sigma = {f \over \sqrt{R_i}}$
Calculated in the model as $\cancel {\sigma = {max(|f|,\sqrt{c_r2\beta}) \over \sqrt{R_i}}}$ Just using $\sigma = {f \over \sqrt{R_i}}$
with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.
##### Plot max(sqrt(c_r*2*beta),|f|)
orig_ecco_def = pyresample.geometry.GridDefinition(lons=lon_orig, lats=lat_orig)
####grate_nearest = pyresample.kd_tree.resample_nearest(orig_def, grate, \
##### targ_def, radius_of_influence=500000, fill_value=None)
####c_rossby_nearest = pyresample.kd_tree.resample_nearest(orig_def, c_rossby, \
#### targ_def, radius_of_influence=500000, fill_value=None)
####btp_nearest = pyresample.kd_tree.resample_nearest(orig_def, btp, \
#### targ_def, radius_of_influence=500000, fill_value=None)
####fcor_nearest = pyresample.kd_tree.resample_nearest(orig_def, fcor, \
#### targ_def, radius_of_influence=500000, fill_value=None)
####tinv_nearest = pyresample.kd_tree.resample_nearest(orig_def, tinv, \
#### targ_def, radius_of_influence=500000, fill_value=None)
####
####tinv_calc=np.sqrt(2*btp_nearest*c_rossby_nearest)
####maxf=np.maximum(tinv_calc,np.abs(fcor_nearest))
absf=np.abs(fcor_nearest)
fig3=MapContourf(lon_targcyc,lat_targcyc,absf,
figsize=(14,5),
# title=" maxf = $max(|f|,\sqrt{c_r*2*B})$")
title="$|f|$")
##### Plot sqrt(R_i)
grate_nearest = pyresample.kd_tree.resample_nearest(orig_def, grate, \
targ_def, radius_of_influence=500000, fill_value=None)
print grate_nearest.min(),grate_nearest.max()
#clevs=np.linspace(0,2200,23)
clevs=np.logspace(0,6,20)
norm=LogNorm()
sqrtgrate_nearest=np.sqrt(grate_nearest)
fig4=MapContourf(lon_targcyc,lat_targcyc,sqrtgrate_nearest,
levels=clevs,
norm=norm,
figsize=(14,5),
title="$\sqrt{R_i}$ (log scale)")
##### Plot max(sqrt(c_r*2*beta),|f|)/sqrt(R_i)
sigma_avg_calc_day=(absf/sqrtgrate_nearest)*86400.
clevs=np.logspace(-3,0,30)
norm=LogNorm()
fig1=MapContourf(lon_targcyc,lat_targcyc,sigma_avg_calc_day,
levels=clevs,
norm=norm,
figsize=(14,5),
extend='both',
title="Eady inverse time scale calc $(days^{-1}) \sim {|f| \over \sqrt{R_i}}$ (log scale)")
clevs=arange(.005,.0625,.0025)
fig=MapContourf(lon_targcyc,lat_targcyc,sigma_avg_calc_day,
levels=clevs,
figsize=(14,5),
extend='both',
title="Eady inverse time scale calc $(days^{-1}) \sim {|f| \over \sqrt{R_i}}$")
maxf=np.maximum(np.sqrt(c_rossby_nearest*2.*btp_nearest),absf)
sigma_maxf_calc_day=(maxf/sqrtgrate_nearest)*86400.
clevs=arange(.005,.0625,.0025)
#fig=MapContourf(lon_targcyc,lat_targcyc,sigma_maxf_calc_day,
# levels=clevs,
# figsize=(14,5),
# extend='both',
# title="Eady inverse time scale calc maxf $(days^{-1}) \sim {|f| \over \sqrt{R_i}}$")
tulloch_eady=Image(filename=basedir+'/steering/tulloch_eady.png')
display(tulloch_eady)
Alternate Sigma Calculation
$ R_i = {f^2 N^2 \over \underbrace{{ {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial y})^2 + {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial z})^2} }_\text{m^4}}\quad\quad = \quad\quad {f^2N^2 \over m^4} $
$\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}$
clevs=arange(.005,.0625,.0025)
fig=MapContourf(lon_targcyc,lat_targcyc,sigma_avg_calc_day,
levels=clevs,
figsize=(14,5),
extend='both',
title="Old Eady inverse time scale calc $(days^{-1}) \sim {|f| \over \sqrt{R_i}}$")
ncalt = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nmonth1.0249-02.nc')
#print(nc.variables)
sigma_avg1 = ncalt.variables['SIGMA_AVG'][0,:,:]
sigma_avg1_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_avg1, \
targ_def, radius_of_influence=500000, fill_value=None)
sigma_avg1_monthavg_day=sigma_avg1_nearest*86400.
clevs=arange(.005,.0625,.0025)
fig=MapContourf(lon_targcyc,lat_targcyc,sigma_avg1_monthavg_day,
levels=clevs,
figsize=(14,5),
extend='both',
title="New Eady inverse time scale calc using alt sigma $(days^{-1}) \sim {{m^2} \over N}$")
bates_len_vel=Image(filename=basedir+'/steering/bates_len_vel.jpg')
display(bates_len_vel)
Zonal Eddy phase speed:
$\cancel{c = - \beta * {L_r^2}}$ $L_r$ too high at equator
$c = - \beta * L_{eddy}^2$
#beta
clevs=arange(1.e-14,25.e-14,1e-14)
norm=matplotlib.colors.Normalize(),
fig12=MapContourf(lon_targcyc,lat_targcyc,btp_nearest,
levels=clevs,
# norm=norm,
addzonal=True,
figsize=(14,5),
extend='both',
title="Beta (1/cm*s)")
clevs=np.logspace(1,15,30)
norm=LogNorm()
l_eddy_calc_sq=l_eddy_calc*l_eddy_calc
fig8=MapContourf(lon_targcyc,lat_targcyc,l_eddy_calc_sq,
addzonal=True,
levels=clevs,
norm=norm,
figsize=(14,5),
extend='both',
title="$L_{eddy}^2 (cm^2)$ log scale")
c=btp_nearest*l_eddy_calc*l_eddy_calc # units=1/(cm s) * cm^2 = cm/s
#clevs=np.linspace(-10,30,10)
clevs=arange(-10,30,1)
norm=LogNorm()
fig12=MapContourf(lon_targcyc,lat_targcyc,c,
addzonal=True,
levels=clevs,
# norm=norm,
figsize=(14,5),
contourlines=True,
extend='both',
title="eddy phase speed c = beta*leddy**2 (cm/s) linear scale(NOTE I'm plotting positive c here - changed to negative after this plot")
c=-1.*c
Hughes Phase Speed (cm/s) from Tulloch Marshall Smith '09
hughes_c=Image(filename=basedir+'/steering/Hughes_phase_speed.png')
display(hughes_c)
Eddy Velocity:
$u_{rms} = alpha*\sigma*L_{eddy}$
############################
##### PLOT L_Eddy
############################
cm2m=.01
alpha=5.
sigma_calc=(absf/sqrtgrate_nearest)
#sigma_calc=(maxf/sqrtgrate_nearest)
u_rms_leddy = alpha*sigma_calc*l_eddy_calc # units = 1*1/s*cm = cm/s
u_rms_leddy_mps=u_rms_leddy*cm2m
clevs=np.logspace(0,6,20)
cm2m=.01
l_eddy_calc_m=l_eddy_calc*cm2m
norm=LogNorm()
fig9=MapContourf(lon_targcyc,lat_targcyc,l_eddy_calc_m,
addzonal=True,
levels=clevs,
norm=norm,
figsize=(14,5),
# extend='both',
title="$L_{eddy} (m)$ log scale")
pylab.xlim([0,3e5])
################# PLOT u_rms with Leddy #######################
u_rms_leddy_mps=u_rms_leddy*.01
clevs=np.linspace(0,.4,40)
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_mps,
addzonal=True,
levels=clevs,
figsize=(14,5),
extend='both',
title="$u_{rms} ({m \over s})$ (alpha=%i) %i * old $\sigma*l_{eddy}$"%(alpha,alpha))
pylab.xlim([0,.2])
sigma_calc_new=sigma_avg1_nearest
u_rms_leddy_new = alpha*sigma_calc_new*l_eddy_calc # units = 1*1/s*cm = cm/s
u_rms_leddy_new_mps=u_rms_leddy_new*cm2m
clevs=np.linspace(0,.4,40)
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_new_mps,
addzonal=True,
levels=clevs,
figsize=(14,5),
extend='both',
title="$u_{rms} ({m \over s})$ (alpha=%i) %i * alt $\sigma*l_{eddy}$"%(alpha,alpha))
pylab.xlim([0,.2])
bates_urms=Image(filename=basedir+'/steering/bates2013-urms.png')
display(bates_urms)
################# PLOT u_rms^2 convert m^2/s^2 #######################
u_rms_leddy_mps_sq=u_rms_leddy_mps*u_rms_leddy_mps
#clevs=np.linspace(0,.06,50)
clevs=arange(0,.06,.001)
#clevs=np.logspace(-5,1,30)
norm=LogNorm()
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_mps_sq,
addzonal=True,
levels=clevs,
# norm=norm,
figsize=(14,5),
extend='both',
title="$u_{rms}^2 ({m^2 \over s^2})$ (alpha=%i) %i * old $\sigma*l_{eddy}$ "%(alpha,alpha))
pylab.xlim([0,.06])
u_rms_leddy_new_mps_sq=u_rms_leddy_new_mps*u_rms_leddy_new_mps
#clevs=np.linspace(0,.06,50)
clevs=arange(0,.06,.001)
#clevs=np.logspace(-5,1,30)
norm=LogNorm()
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_new_mps_sq,
addzonal=True,
levels=clevs,
# norm=norm,
figsize=(14,5),
extend='both',
title="$u_{rms}^2 ({m^2 \over s^2})$ (alpha=%i) %i * Alt $\sigma*l_{eddy}$ "%(alpha,alpha))
pylab.xlim([0,.06])
urms_bates=Image(filename=basedir+'/steering/Bates_urms_sq.jpg')
display(urms_bates)
################# PLOT (u-c)^2 convert m^2/s^2 #######################
umean10_nearest = pyresample.kd_tree.resample_nearest(orig_def, umean10,
targ_def, radius_of_influence=500000, fill_value=None)
umean10_mps=umean10_nearest*1.e-2
clevs=np.linspace(-60,60,40)
fig11=MapContourf(lon_targcyc,lat_targcyc,umean10_nearest,
addzonal=True,
figsize=(14,5),
levels=clevs,
extend='both',
title="U zonal velocity (cm/s)")
pylab.xlim([-20.,20.])
# get 2D versions of the lat and lon variables add longitude start point here!
lon_ecco_orig, lat_ecco_orig = np.meshgrid(lon_ecco[:], lat_ecco[:])
orig_ecco_def = pyresample.geometry.GridDefinition(lons=lon_ecco_orig, lats=lat_ecco_orig)
ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_ecco_def, ecco_u, \
targ_def, radius_of_influence=500000, fill_value=None)
ecco_u_nearest_masked=np.ma.masked_where(ecco_u_nearest<-10, ecco_u_nearest, copy=True)
#print ecco_u_cyc_masked[:,0]
clevs=np.linspace(-60,60,40)
#ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_def, ecco_u, \
# targ_def, radius_of_influence=500000, fill_value=None)
fig11=MapContourf(lon_targcyc,lat_targcyc,ecco_u_nearest_masked*100,
addzonal=True,
figsize=(14,5),
levels=clevs,
# extend='both',
title="ECCO U zonal velocity (cm/s)")
pylab.xlim([-20.,20.])
clevs=-np.logspace(3,-3,20)
print clevs
norm=matplotlib.colors.SymLogNorm(linthresh=1e-4, vmin=-1000, vmax=-.001)
fig12=MapContourf(lon_targcyc,lat_targcyc,c,
addzonal=True,
levels=clevs,
norm=norm,
figsize=(14,5),
# extend='both',
title="eddy phase speed c = beta*leddy**2 (cm/s) log scale")
pylab.xlim([-200.,0.])
uminusc=umean10_nearest-c #units cm/s
#clevs=np.logspace(-2,3,10)
clevs=np.linspace(0,25,25)
#norm=LogNorm()
fig13=MapContourf(lon_targcyc,lat_targcyc,uminusc,
addzonal=True,
# norm=norm,
levels=clevs,
figsize=(14,5),
# extend='both',
title="U-c (cm/s) (max = %i/min = %i) (c calculated using Leddy)"%(uminusc.max(),uminusc.min()))
pylab.xlim([0,25])
bates_uminusc=Image(filename=basedir+'/steering/bates2013_uminusc.png')
display(bates_uminusc)
uminusc_mps_sq=uminusc*uminusc*1.e-4 #units cm/s*cm/s*m/100cm*m/100cm
clevs=np.linspace(0,.06,41)
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_mps_sq,
addzonal=True,
levels=clevs,
figsize=(14,5),
# extend='both',
title=" ${(U-c)}^2 {(m^2/s^2)}$ max = %i/min = %i (c calculated using Leddy) "%(uminusc_mps_sq.max(),uminusc_mps_sq.min()))
bates_uminuscsq=Image(filename=basedir+'/steering/Bates_u_minus_c_sq.jpg')
display(bates_uminuscsq)
uminusc_sq_ovr_urms_sq = uminusc_mps_sq/u_rms_leddy_mps_sq #units cm/s / cm/s
uminusc_sq_ovr_urms_sq_masked=np.ma.masked_where(uminusc_sq_ovr_urms_sq >30., uminusc_sq_ovr_urms_sq, copy=True)
clevs=arange(0,5,.1)
#clevs=np.logspace(-3,3,21)
norm=LogNorm()
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_sq,
addzonal=True,
levels=clevs,
# norm=norm,
figsize=(14,5),
extend='both',
title="${(U-c)}^2 \over u_{rms}^2$")
pylab.xlim([0,8])
uminusc_sq_ovr_urms_new_sq = uminusc_mps_sq/u_rms_leddy_new_mps_sq #units cm/s / cm/s
clevs=arange(0,5,.1)
#clevs=np.logspace(-3,3,21)
norm=LogNorm()
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_new_sq,
addzonal=True,
levels=clevs,
# norm=norm,
figsize=(14,5),
extend='both',
title="${(U-c)}^2 \over u_{rms}^2$ using Alt sigma")
pylab.xlim([0,8])
clevs=np.logspace(-3,4,50)
norm=LogNorm()
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_new_sq,
addzonal=True,
levels=clevs,
norm=norm,
figsize=(14,5),
extend='both',
title="${(U-c)}^2 \over u_{rms}^2$ max=%i/min=%i (using Alt sigma) log scale"%(uminusc_sq_ovr_urms_new_sq.max(),uminusc_sq_ovr_urms_new_sq.min()))
pylab.xlim([0,100])
bates_uminusc_sq_ovr_urms_sq =Image(filename=basedir+'/steering/Bates_suppress.jpg')
display(bates_uminusc_sq_ovr_urms_sq )
$K=u_{rms}∗L_{mix}$
$L_{mix} = \Gamma * L_{eddy} * Suppression$
$Suppression= {1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2 )}$
$L_{mix} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}$
cm2m=.01
u_rms_leddy_mps=u_rms_leddy*cm2m
u_rms_leddy_mps_sq=u_rms_leddy_mps*u_rms_leddy_mps
gamma=.35
b1=4.0
supp=1./(1+b1*uminusc_sq_ovr_urms_sq)
clevs=np.linspace(0,1,30)
clevs=arange(0,1,.05)
norm=matplotlib.colors.Normalize()
fig=MapContourf(lon_targcyc,lat_targcyc,supp,
addzonal=True,
norm=norm,
levels=clevs,
figsize=(14,5),
title="Suppression factor ${1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2)}$ using old sigma")
pylab.xlim([0,.8])
supp_new=1./(1+b1*uminusc_sq_ovr_urms_new_sq)
clevs=np.linspace(0,1,20)
norm=matplotlib.colors.Normalize()
fig=MapContourf(lon_targcyc,lat_targcyc,supp_new,
addzonal=True,
norm=norm,
levels=clevs,
figsize=(14,5),
title="Suppression factor ${1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2)}$ using Alt sigma")
pylab.xlim([0,.8])
bates_suppression =Image(filename=basedir+'/steering/bates2013_suppression.png')
display(bates_suppression )
gammax5=gamma*5
lmix_new_m=gammax5*l_eddy_calc_m*supp_new
clevs=np.logspace(0,6,20)
norm=LogNorm()
fig9=MapContourf(lon_targcyc,lat_targcyc,l_eddy_calc_m,
addzonal=True,
levels=clevs,
norm=norm,
figsize=(14,5),
# extend='both',
title="$L_{eddy} (m)$ log scale")
pylab.xlim([0,3e5])
clevs=np.linspace(0,100e3,10)
norm=matplotlib.colors.Normalize()
fig=MapContourf(lon_targcyc,lat_targcyc,lmix_new_m,
addzonal=True,
norm=norm,
levels=clevs,
figsize=(14,5),
title="lmix (m) = gamma*leddy*suppression (using Alt sigma ... also used gamma x 5 (5*.35))")
#pylab.xlim([0,10e3])
clevs=np.linspace(0,.4,40)
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_new_mps,
addzonal=True,
levels=clevs,
figsize=(14,5),
extend='both',
title="$u_{rms} ({m \over s})$ (alpha=%i) %i * alt $\sigma*l_{eddy}$"%(alpha,alpha))
pylab.xlim([0,.2])
k_new=u_rms_leddy_new_mps*lmix_new_m
clevs=np.linspace(0,10e3,40)
norm=matplotlib.colors.Normalize()
fig=MapContourf(lon_targcyc,lat_targcyc,k_new,
addzonal=True,
norm=norm,
levels=clevs,
figsize=(14,5),
title="k=urms*lmix (using Alt sigma for calculation of urms and gamma * 5 for calculation of lmix")
pylab.xlim([0,10e3])
urms_bates2=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-002.jpg')
urms_bates3=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-003.jpg')
urms_bates4=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-004.jpg')
urms_bates5=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-005.jpg')
urms_bates6=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-006.jpg')
urms_bates7=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-007.jpg')
urms_bates8=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-008.jpg')
urms_bates9=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-009.jpg')
urms_bates10=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-010.jpg')
urms_bates11=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-011.jpg')
display(urms_bates3,urms_bates4,urms_bates5,urms_bates6,urms_bates7,urms_bates8,urms_bates9,urms_bates10,urms_bates11)
CESM $U_{mean} = U_{resolved}$ integraged over {10,50,100,200,500} meters
umean_cesm=Image(filename=basedir+'/steering/umean.jpg')
display(umean_cesm)
sigma_avg_calc=(absf/sqrtgrate_nearest)
rhines_km=sigma_avg_calc/btp_nearest*cm2km
print rhines_km.max(),rhines_km.min()
clevs=[10,20,30,40,50,60,80,100,150,225]
clevs=arange(0,225,10)
fig=MapContourf(lon_targcyc,lat_targcyc,rhines_km,
addzonal=True,
levels=clevs,
figsize=(14,5),
extend='both',
title="Rhines Scale $(km) \sim {Sigma \over B}$")
popdata = Dataset(tfile, 'r')
print tfile
print popdata
shear_max = popdata.variables('SHEAR_MAX')
from ctd import plot_section
from oceans.colormaps import cm
#from oceans.datasets import etopo_subset
def laplace_X(F, M):
"""1D Laplace Filter in X-direction.
http://www.trondkristiansen.com/wp-content/uploads/downloads/
2010/09/laplaceFilter.py"""
jmax, imax = F.shape
# Add strips of land.
F2 = np.zeros((jmax, imax + 2), dtype=F.dtype)
F2[:, 1:-1] = F
M2 = np.zeros((jmax, imax + 2), dtype=M.dtype)
M2[:, 1:-1] = M
MS = M2[:, 2:] + M2[:, :-2]
FS = F2[:, 2:] * M2[:, 2:] + F2[:, :-2] * M2[:, :-2]
return np.where(M > 0.5, (1 - 0.25 * MS) * F + 0.25 * FS, F)
def laplace_Y(F, M):
"""1D Laplace Filter in Y-direction.
http://www.trondkristiansen.com/wp-content/uploads/downloads/
2010/09/laplaceFilter.py"""
jmax, imax = F.shape
# Add strips of land.
F2 = np.zeros((jmax + 2, imax), dtype=F.dtype)
F2[1:-1, :] = F
M2 = np.zeros((jmax + 2, imax), dtype=M.dtype)
M2[1:-1, :] = M
MS = M2[2:, :] + M2[:-2, :]
FS = F2[2:, :] * M2[2:, :] + F2[:-2, :] * M2[:-2, :]
return np.where(M > 0.5, (1 - 0.25 * MS) * F + 0.25 * FS, F)
def laplace_filter(F, M=None):
"""Laplace filter a 2D field with mask. The mask may cause laplace_X and
laplace_Y to not commute. Take average of both directions.
http://www.trondkristiansen.com/wp-content/uploads/downloads/
2010/09/laplaceFilter.py"""
if not M:
M = np.ones_like(F)
return 0.5 * (laplace_X(laplace_Y(F, M), M) +
laplace_Y(laplace_X(F, M), M))
def pop_subset(var='',llcrnrlon=None, urcrnrlon=None, llcrnrlat=None, urcrnrlat=None, tfile='', smoo=False, subsample=False):
### Read Pop file
popdata = Dataset(tfile, 'r')
popvar = popdata.variables[var][0,:,:,:]
print popvar.shape
popvar=np.rollaxis(popvar, 1)
print popvar.shape
popvar=np.rollaxis(popvar, 2,1)
print popvar.shape
### Get lats lons for original and target grids for remapping
### Process input lons for pyregrid (-180:180) and output lons for matplotlib idiosyncracies
lon_orig = popdata.variables["TLONG"][:]
lat_orig = popdata.variables["TLAT"][:]
lev_targ = popdata.variables["z_t"][:]
# setup for using pyresample to interpolate plot
# make sure lon is in range of -180:180 for pyregrid
lon_orig = (lon_orig + 180) % (360) - 180
# open dataset for output grid get latout,lonout
fc = Dataset(basedir+'/steering/co2flux_fossil_1751-2006-monthly_0.9x1.25_c20100204.nc')
lon= fc.variables['lon'][:]
lat = fc.variables['lat'][:]
fc.close()
# get 2D versions of the lat and lon variables add longitude start point here!
lon_targ, lat_targ = np.meshgrid(lon, -1.*lat)
# add cyclic point to output lat,lon arrays for plot routines
lon_targcyc = np.concatenate((lon_targ, lon_targ[:, 0:1]), axis=1)
lat_targcyc = np.concatenate((lat_targ, lat_targ[:, 0:1]), axis=1)
# make sure lon is in range of -180:180 for pyregrid
lon_targcyc = (lon_targcyc + 180) % (360) - 180
#make longitudes monotonic for contouring add longitude period to end
lon_targ = np.where(np.greater_equal(lon_targcyc,min(lon_targcyc[:,0])),lon_targcyc-360,lon_targcyc)+360.
endlon=lon_targ.shape[1]-1
lon_targ[:,endlon]=lon_targ[:,endlon]+360.
lat_targ=lat_targcyc
#When all the necessary inputs are read in memory, that's where we start using pyresample's magic.
#Create a pyresample object holding the origin (POP) grid:
orig_def = pyresample.geometry.GridDefinition(lons=lon_orig, lats=lat_orig)
#Create another pyresample object for the target (curvilinear) grid:
targ_def = pyresample.geometry.GridDefinition(lons=lon_targcyc, lats=lat_targcyc)
print lon_orig.shape,lat_orig.shape,lon_targ.shape,lat_targ.shape,popvar.shape
popvar_remap = pyresample.kd_tree.resample_nearest(orig_def, popvar, \
targ_def, radius_of_influence=500000, fill_value=None)
# Select data subset.
maskx = np.logical_and(lon_targ >= llcrnrlon, lon_targ <= urcrnrlon)
masky = np.logical_and(lat_targ >= llcrnrlat, lat_targ <= urcrnrlat)
maskz = np.array([z in depths for z in depth])
### set region limits
###res = get_indices(llcrnrlat, urcrnrlat, llcrnrlon, urcrnrlon, lon_targ, lat_targ)
###lon, lat = np.meshgrid(lons[res[0]:res[1]], lats[res[2]:res[3]])
####popvar_clip = ncdataset.variables[var][int(res[2]):int(res[3]),int(res[0]):int(res[1])]
# if smoo:
# popvar_remap = laplace_filter(popvar_remap, M=None)
if subsample:
lon_targ, lat_targ, popvar_remap = lon_targ[::subsample], lat_targ[::subsample], popvar_remap[::subsample]
return lon_targ, lat_targ, lev_targ, popvar_remap
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from pandas import Panel4D, Panel
from netCDF4 import Dataset, num2date
from ..ff_tools import get_profile, wrap_lon180, wrap_lon360
#
def map_limits(m):
lons, lats = wrap_lon360(m.boundarylons), m.boundarylats
boundary = dict(llcrnrlon=min(lons),
urcrnrlon=max(lons),
llcrnrlat=min(lats),
urcrnrlat=max(lats))
return boundary
longitude = -25.5
llcrnrlon=longitude
urcrnrlon=longitude
llcrnrlon, urcrnrlon = map(wrap_lon360, (llcrnrlon, urcrnrlon))
lon = wrap_lon360(nc.variables.pop('lon')[:])
lat = nc.variables.pop('lat')[:]
depth = nc.variables.pop('depth')[:]
times = nc.variables.pop('time')
times = num2date(times[:], times.units, calendar='365_day')
times = [time.strftime('%b') for time in times]
longitude = -25.5
llcrnlon=longitude
tfile=basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.0249-01-01.nc'
nc= Dataset(tfile, 'r')
lon = wrap_lon360(nc.variables.pop('lon')[:])
print lon
# Select data subset.
#maskx = np.logical_and(lon >= llcrnrlon, lon <= urcrnrlon)
#masky = np.logical_and(lat >= llcrnrlat, lat <= urcrnrlat)
#maskz = np.array([z in depths for z in depth])
levels=slice(0,20)
depths = [0, 10, 20, 30, 50, 75, 100, 125, 150, 200, 250, 300, 400, 500,
600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1750,
2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, 7000,
7500, 8000, 8500, 9000][levels]
print depths
var='temperature'
v = dict(temperature='t', dissolved_oxygen='o', salinity='s',
oxygen_saturation='O', apparent_oxygen_utilization='A',
phosphate='p', silicate='p', nitrate='n')
start = '%s_' % v[var]
print start
variables = dict()
for variable in nc.variables.keys():
if variable.startswith(start):
subset = nc.variables[variable][..., maskz, masky, maskx]
data = Panel4D(subset, major_axis=lat, minor_axis=lon,
labels=np.atleast_1d(times),
items=np.atleast_1d(depth))
variables.update({d[variable]: data})
print variables
#dataset = dataset.clip(27.5, 37.5)
#levels = np.arange(dataset.min().min(), dataset.max().max(), 0.02)
#dataset.lat = dataset.columns.values.astype(float)
#dataset.lon = len(dataset.lat) * [longitude]
tfile=basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.0249-01-01.nc'
longitude=-25.5
lons, lats, levs, shear = pop_subset(var='SHEAR_MAX',
llcrnrlon=longitude,
urcrnrlon=longitude,
llcrnrlat=-70,
urcrnrlat=70,
smoo=True,
tfile=tfile)
print shear.shape
print shear[:,240,:].shape
print levs.shape
fig, ax = plt.subplots()
latm, levm = np.meshgrid(lats[:,240], levs[::-1])
cs = ax.pcolormesh(latm, levm, shear[:,240,:].transpose())
fig, ax, cbar = plot_section(shear, cmap=ocm.avhrr, marker=None,
levels=levels, figsize=(10, 4))
X = ax.get_xticks()
offset = 1.01
ax.set_xlim(X.min() * offset, X.max() * offset)
new_labels = np.linspace(dataset.lat.min(), dataset.lat.max(), len(X))
new_labels = [u'%i\u00B0' % num for num in new_labels]
ax.set_xticklabels(new_labels)
ax.set_xlabel(u'WOA09 Salinity Section at Longitude %3.1f\u00B0' % longitude)
ax.set_ylabel('Depth [m]')
axin = inset_axes(ax, width="40%", height="40%", loc=4)
inmap = Basemap(projection='ortho', lon_0=longitude, lat_0=0,
ax=axin, anchor='NE')
inmap.bluemarble()
_ = inmap.plot(dataset.lon, dataset.lat, 'r', latlon=True)